home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpkdrv.zip / SS.FOR < prev    next >
Text File  |  1984-01-07  |  19KB  |  656 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL SSITS(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SSITS(LUNIT)
  14. C     LUNIT IS THE OUTPUT UNIT NUMBER
  15. C
  16. C     TESTS
  17. C        SSICO,SSIFA,SSISL,SSIDI,SSPCO,SSPFA,SSPSL,SSPDI
  18. C
  19. C     LINPACK. THIS VERSION DATED 08/14/78 .
  20. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  21. C
  22. C     SUBROUTINES AND FUNCTIONS
  23. C
  24. C     LINPACK SSICO,SSISL,SSIDI,SSPCO,SSPSL,SSPDI
  25. C     EXTERNAL SSIXX,SMACH
  26. C     BLAS SAXPY,SSWAP,SASUM,SCOPY
  27. C     FORTRAN ABS,AMAX1,FLOAT,IABS
  28. C
  29. C     INTERNAL VARIABLES
  30. C
  31.       REAL A(15,15),AINV(15,15),ASAVE(15,15),AJK,AJKP1
  32.       REAL AP(120),B(15),X(15),XEXACT(15)
  33.       REAL XP(15),T,Z(15)
  34.       REAL ANORM,AINORM,COND,COND1,DET(2),DETP(2)
  35.       REAL EN,ENORM,EPS,FNORM,ONEPX,Q(6),QS(6),RCOND
  36.       REAL RCONDP,RNORM,S,SASUM,SMACH,XNORM
  37.       INTEGER I,INERT(3),INERTP(3),IQ(6),J,JB
  38.       INTEGER K,KASE,KSING,KM1,KOUNT,KPFAIL,KPVT(15),KPVTP(15)
  39.       INTEGER KS,KSTEP,KSUSP(6),LDA,LUNIT,N,NPRINT
  40.       LOGICAL KPF
  41. C
  42.       LDA = 15
  43. C
  44. C     WRITE MATRIX AND SOLUTIONS IF  N .LE. NPRINT
  45. C
  46.       NPRINT = 3
  47. C
  48.       WRITE (LUNIT,540)
  49.       WRITE (LUNIT,950)
  50. C
  51.       DO 10 I = 1, 6
  52.          KSUSP(I) = 0
  53.    10 CONTINUE
  54.       KSING = 0
  55.       KPFAIL = 0
  56. C
  57. C     SET EPS TO ROUNDING UNIT FOR REAL ARITHMETIC
  58. C
  59.       EPS = SMACH(1)
  60.       WRITE (LUNIT,550) EPS
  61.       WRITE (LUNIT,530)
  62. C
  63. C     START MAIN LOOP
  64. C
  65.       KASE = 1
  66.    20 CONTINUE
  67. C
  68. C        GENERATE TEST MATRIX
  69. C
  70.          CALL SSIXX(A,LDA,N,KASE,LUNIT)
  71. C
  72. C        N = 0 SIGNALS NO MORE TEST MATRICES
  73. C
  74. C     ...EXIT
  75.          IF (N .LE. 0) GO TO 520
  76.          ANORM = 0.0E0
  77.          DO 30 J = 1, N
  78.             ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
  79.    30    CONTINUE
  80.          WRITE (LUNIT,700) ANORM
  81. C
  82.          IF (N .GT. NPRINT) GO TO 50
  83.             WRITE (LUNIT,530)
  84.             DO 40 I = 1, N
  85.                WRITE (LUNIT,740) (A(I,J), J = 1, N)
  86.    40       CONTINUE
  87.             WRITE (LUNIT,530)
  88.    50    CONTINUE
  89. C
  90. C        GENERATE EXACT SOLUTION
  91. C
  92.          XEXACT(1) = 1.0E0
  93.          IF (N .GE. 2) XEXACT(2) = 0.0E0
  94.          IF (N .LE. 2) GO TO 70
  95.             DO 60 I = 3, N
  96.                XEXACT(I) = -XEXACT(I-2)
  97.    60       CONTINUE
  98.    70    CONTINUE
  99. C
  100. C        SAVE MATRIX AND GENERATE R.H.S.
  101. C
  102.          DO 90 I = 1, N
  103.             B(I) = 0.0E0
  104.             DO 80 J = 1, N
  105.                ASAVE(I,J) = A(I,J)
  106.                B(I) = B(I) + A(I,J)*XEXACT(J)
  107.    80       CONTINUE
  108.             X(I) = B(I)
  109.             XP(I) = X(I)
  110.    90    CONTINUE
  111. C
  112. C        FACTOR AND ESTIMATE CONDITION
  113. C
  114.          RCOND = -1.0E0
  115.          CALL SSICO(A,LDA,N,KPVT,RCOND,Z)
  116.          WRITE (LUNIT,820) (KPVT(I), I = 1, N)
  117. C
  118. C        OUTPUT NULL VECTOR
  119. C
  120.          IF (N .GT. NPRINT) GO TO 110
  121.             WRITE (LUNIT,760)
  122.             DO 100 I = 1, N
  123.                WRITE (LUNIT,770) Z(I)
  124.   100       CONTINUE
  125.             WRITE (LUNIT,530)
  126.   110    CONTINUE
  127. C
  128. C        FACTOR PACKED FORM AND COMPARE
  129. C
  130.          KPF = .FALSE.
  131.          K = 0
  132.          DO 130 J = 1, N
  133.             DO 120 I = 1, J
  134.                K = K + 1
  135.                AP(K) = ASAVE(I,J)
  136.   120       CONTINUE
  137.   130    CONTINUE
  138.          RCONDP = -1.0E0
  139.          CALL SSPCO(AP,N,KPVTP,RCONDP,Z)
  140.          IF (RCONDP .EQ. RCOND) GO TO 140
  141.             WRITE (LUNIT,840)
  142.             WRITE (LUNIT,880) RCOND,RCONDP
  143.             KPF = .TRUE.
  144.   140    CONTINUE
  145.          K = 0
  146.          KOUNT = 0
  147.          DO 160 J = 1, N
  148.             DO 150 I = 1, J
  149.                K = K + 1
  150.                IF (AP(K) .NE. A(I,J)) KOUNT = KOUNT + 1
  151.   150       CONTINUE
  152.   160    CONTINUE
  153.          IF (KOUNT .EQ. 0) GO TO 170
  154.             WRITE (LUNIT,840)
  155.             WRITE (LUNIT,890) KOUNT
  156.             KPF = .TRUE.
  157.   170    CONTINUE
  158. C
  159. C        TEST FOR SINGULARITY
  160. C
  161.          IF (RCOND .GT. 0.0E0) GO TO 180
  162.             WRITE (LUNIT,750) RCOND
  163.             KSING = KSING + 1
  164.             CALL SSIDI(A,LDA,N,KPVT,DET,INERT,Z,110)
  165.             WRITE (LUNIT,790) DET(1)
  166.             WRITE (LUNIT,810) INERT
  167.          GO TO 510
  168.   180    CONTINUE
  169.             COND = 1.0E0/RCOND
  170.             WRITE (LUNIT,570) COND
  171.             ONEPX = 1.0E0 + RCOND
  172.             IF (ONEPX .EQ. 1.0E0) WRITE (LUNIT,560)
  173. C
  174. C           COMPUTE INVERSE, DETERMINANT, INERTIA, TRUE CONDITION
  175. C
  176.             DO 200 J = 1, N
  177.                DO 190 I = 1, J
  178.                   AINV(I,J) = A(I,J)
  179.   190          CONTINUE
  180.   200       CONTINUE
  181.             CALL SSIDI(AINV,LDA,N,KPVT,DET,INERT,Z,111)
  182.             AINORM = 0.0E0
  183.             DO 220 J = 1, N
  184.                DO 210 I = J, N
  185.                   AINV(I,J) = AINV(J,I)
  186.   210          CONTINUE
  187.                AINORM = AMAX1(AINORM,SASUM(N,AINV(1,J),1))
  188.   220       CONTINUE
  189.             COND1 = ANORM*AINORM
  190.             WRITE (LUNIT,580) COND1
  191.             WRITE (LUNIT,790) DET(1)
  192.             WRITE (LUNIT,800) DET(2)
  193.             WRITE (LUNIT,810) INERT
  194. C
  195. C           SOLVE A*X = B
  196. C
  197.             CALL SSISL(A,LDA,N,KPVT,X)
  198. C
  199.             IF (N .GT. NPRINT) GO TO 240
  200.                WRITE (LUNIT,590)
  201.                DO 230 I = 1, N
  202.                   WRITE (LUNIT,780) X(I)
  203.   230          CONTINUE
  204.                WRITE (LUNIT,530)
  205.   240       CONTINUE
  206. C
  207. C           MORE PACKED COMPARE
  208. C
  209.             CALL SSPSL(AP,N,KPVTP,XP)
  210.             KOUNT = 0
  211.             DO 250 I = 1, N
  212.                IF (XP(I) .NE. X(I)) KOUNT = KOUNT + 1
  213.   250       CONTINUE
  214.             IF (KOUNT .EQ. 0) GO TO 260
  215.                WRITE (LUNIT,840)
  216.                WRITE (LUNIT,900) KOUNT
  217.                KPF = .TRUE.
  218.   260       CONTINUE
  219.             CALL SSPDI(AP,N,KPVTP,DETP,INERTP,Z,111)
  220.             IF (DETP(1) .EQ. DET(1) .AND. DETP(2) .EQ. DET(2))
  221.      *         GO TO 270
  222.                WRITE (LUNIT,840)
  223.                WRITE (LUNIT,910) DETP
  224.                KPF = .TRUE.
  225.   270       CONTINUE
  226.             KOUNT = 0
  227.             K = 0
  228.             DO 290 J = 1, N
  229.                DO 280 I = 1, J
  230.                   K = K + 1
  231.                   IF (AP(K) .NE. AINV(I,J)) KOUNT = KOUNT + 1
  232.   280          CONTINUE
  233.   290       CONTINUE
  234.             IF (KOUNT .EQ. 0) GO TO 300
  235.                WRITE (LUNIT,840)
  236.                WRITE (LUNIT,920) KOUNT
  237.                KPF = .TRUE.
  238.   300       CONTINUE
  239. C
  240. C           RECONSTRUCT  A  FROM FACTORS.
  241. C
  242.             K = 1
  243.   310       IF (K .GT. N) GO TO 400
  244.                KM1 = K - 1
  245.                IF (KPVT(K) .LT. 0) GO TO 340
  246. C
  247. C                 1 BY 1
  248. C
  249.                   IF (KM1 .LT. 1) GO TO 330
  250.                   DO 320 JB = 1, KM1
  251.                      J = K - JB
  252.                      AJK = A(K,K)*A(J,K)
  253.                      T = AJK
  254.                      CALL SAXPY(J,T,A(1,K),1,A(1,J),1)
  255.                      A(J,K) = -AJK
  256.   320             CONTINUE
  257.   330             CONTINUE
  258.                   KSTEP = 1
  259.                GO TO 370
  260.   340          CONTINUE
  261. C
  262. C                 2 BY 2
  263. C
  264.                   IF (KM1 .LT. 1) GO TO 360
  265.                   DO 350 JB = 1, KM1
  266.                      J = K - JB
  267.                      AJK = A(K,K)*A(J,K) + A(K,K+1)*A(J,K+1)
  268.                      AJKP1 = A(K,K+1)*A(J,K) + A(K+1,K+1)*A(J,K+1)
  269.                      T = AJK
  270.                      CALL SAXPY(J,T,A(1,K),1,A(1,J),1)
  271.                      T = AJKP1
  272.                      CALL SAXPY(J,T,A(1,K+1),1,A(1,J),1)
  273.                      A(J,K) = -AJK
  274.                      A(J,K+1) = -AJKP1
  275.   350             CONTINUE
  276.   360             CONTINUE
  277.                   KSTEP = 2
  278.   370          CONTINUE
  279. C
  280. C              SWAP
  281. C
  282.                KS = IABS(KPVT(K))
  283.                IF (KS .EQ. K) GO TO 390
  284.                   CALL SSWAP(KS,A(1,KS),1,A(1,K),1)
  285.                   DO 380 JB = KS, K
  286.                      J = K + KS - JB
  287.                      T = A(J,K)
  288.                      A(J,K) = A(KS,J)
  289.                      A(KS,J) = T
  290.   380             CONTINUE
  291.                   IF (KSTEP .EQ. 2)
  292.      *               CALL SSWAP(1,A(KS,K+1),1,A(K,K+